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Abstract 

We discuss artificial boundary conditions for stationary Navier-Stokes flows past bodies in the 
half-plane, for a range of low Reynolds numbers. When truncating the half-plane to a finite domain 
for numerical purposes, artificial boundaries appear. We present an explicit Dirichlet condition for 
the velocity at these boundaries in terms of an asymptotic expansion for the solution to the problem. 
We show a substantial increase in accuracy of the computed values for drag and lift when compared 
with results for traditional boundary conditions. We also analyze the qualitative behavior of the 
solutions in terms of the streamlines of the flow. The new boundary conditions are universal in the 
sense that they depend on a given body only through one constant, which can be determined in a 
feed-back loop as part of the solution process. 
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1 Introduction 

We numerically solve the Navier-Stokes equations for the flow past a body moving at constant veloc- 
ity parallel to the boundary of a half-plane. We are particularly interested in the computation of the 
hydrodynamic forces acting on the body in the case where the body is small, and the flow is laminar. 
We first introduce the viscous length scale, 

iv = — = ~, (1) 

with v the dynamic viscosity of the fluid and Uoo and vb the magnitude of the velocity field at infinity (as 
viewed from the moving body) and the magnitude of the body's translation velocity (as viewed from the 
fluid at rest), respectively. In addition to this dynamic length, there are two geometrical lengths which 
are important in this problem: the body-center-to-wall (hereafter just "body-wall") distance d and the 
body size 2r, d > r. In terms of the viscous length, we define the Reynolds number Re as 

R« = ^ = £. (2) 

We shall keep to small but non-negligible values of the Reynolds number throughout this work, i.e., 
Re = 0.5, . . . , 25. In this regime, it is expected that the flow remains stationary and that viscous and 
inertial phenomena have similar importance, i.e., the flow is neither creeping nor turbulent. 

Since we truncate the unbounded domain to finite sub-domains for the numerical treatment, the ques- 
tion of boundary conditions at the resulting artificial boundaries arises. We show that, when compared to 
traditional methods of "open" boundary conditions (see for example |14|). a significant gain in accuracy 
in the computed values of drag and lift can be obtained when using the asymptotic expansion for the 
velocity field constructed in [2] as Dirichlet boundary conditions. Our method of adaptive boundary con- 
ditions is universal in the sense that it depends on the boundary conditions at the body surface and the 
shape of the body only through one constant. We mainly concentrate on drag and lift for the quantitative 
comparison of different boundary conditions because they are important quantities in engineering and 
theoretical work alike. Besides increased accuracy, the qualitative behavior of the flow within the com- 
putational domain is also improved with our boundary conditions in the sense that the streamlines are 
not significantly influenced by the artificial boundary, contrary to the cases where traditional boundary 
conditions are used. 

We would like to emphasize that the aim of this paper is not to constitute a benchmark for the 
considered problem nor to achieve the highest possible precision, but rather to highlight the fundamental 
importance of the choice of boundary conditions when precision and qualitative correctness of the flow 
patterns are desired alongside a decrease in hardware requirements (especially in memory due to smaller 
computational domains). In order to make the results easily accessible and useful for applications, we 
use for the numerical implementation a widely used commercial code targeted for industrial applications, 
namely COMSOL Multiphasics. See HH1E3] for recent research in CFD using COMSOL Multiphysics in 
low Reynolds numbers regimes. 

The present work is part of an ongoing project in a bottom-up approach to problems with fluid- 
structure interaction, starting with the mathematical analysis of the equations and going all the way 
to numerical applications. Conceptually, we make heavy use of the equivalence between the present 
problem, which we shall refer to as the "body-problem", and a problem without a body, but with a force 
term of compact support, which we shall refer to as the "force-problem". For the mathematical treatment 
of the Oseen- linearized force-problem, see |21j . for an existence theorem for the nonlinear force-problem 
see [22] , and for the proof of uniqueness of solutions and the equivalence of the body-problem and the 
force problem see [23]. A precise bound on the vorticity for the force-problem was derived in [3], which 
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is the key input for the extraction, in [2] , of the asymptotic expansion for the velocity field up to second 
order, which is used in the current work to define the adaptive boundary conditions. For a general 
introduction to the mathematical method used for the derivation of adaptive boundary conditions like 
the ones presented here, see |19| . 

This work is part of an ongoing effort to achieve higher numerical efficiency in the simulation of 
exterior problems. Adaptive boundary conditions of the type described here have been used with success 
for the problem in the full plane, see [4], [5] and (24|) and the three dimensional case is discussed in |18| . 
In three dimensions, other approaches for the numerical treatment of motions near a wall at low Reynolds 
numbers have been put into practice. See for example [7J, where the authors use what we shall refer to 
as classic boundary conditions, and [36J, which involves what we shall call simple boundary conditions. 
Other types of artificial boundary conditions for incompressible viscous flows in exterior domains have 
been developed in the past for various cases. For the time-dependent Oseen-linearized Navier-Stokes flow 
in the full plane, artificial boundary conditions for the normal constraint were proposed in the form of 
differential equations in |16j . For the nonlinear Navier-Stokes equations in the full space, a boundary 
condition based on the leading asymptotic terms of the solution was obtained in [26 . For high Reynolds 
number (Re > 10 7 ) flows in aerodynamics (mainly compressible, but including the incompressible limit) 
artificial boundary conditions based on the Calderon-Ryaben'kii method involving difference potentials 
and pseudo-differential boundary operators has culminated in the work in [32_ , which includes numerical 
applications. Another artificial boundary condition based on a detailed mathematical analysis of bounded 
domains successively approximating an unbounded one is derived in [9] with error estimates given in jS] . A 
heuristic approach to artificial boundary conditions, the so-called "no-boundary" or "do nothing" boundary 
conditions, was pioneered in |251 1271 |2"U] , It consists in extending the governing equations to the artificial 
boundary ([15] showed that this method implicitly imposes the (p + l) th derivative of the solution to 
vanish at a point close to the outflow, where p is the degree of the finite elements), but this approach 
was applied for domains where only the outflow is an artificial boundary, together with prescribed inflow, 
and its validity in the case of full exterior domains is unclear. Finally, in [T7J, a procedure based on the 
conservation of mass and vorticity considerations is used to extrapolate the radial velocity at the artificial 
boundaries for the problem of the flow past a body placed in a channel, with improved behavior of the 
streamlines at the outlet even in non-stationary cases. 

Experimental work on single bodies moving slowly parallel to a wall has been a recurring topic for over 
fifty years, with works such as |13| and [T|. More recently, in [31 and [30J, experiments were reported that 
show that the transverse force on small spheroidal bubbles moving close to a recipient wall appears to 
change sign for particular parameters of the flow, the bubble-fluid interface type and small deformations 
of the bubbles. Two types of bubble-fluid interfaces were studied, so-called "clean" and "contaminated" 
bubbles, which are modeled mathematically as "slip" and "noslip" boundary conditions, respectively. The 
authors suggest that the change in sign of the transverse force is due to two competing mechanisms. On 
the one hand, vorticity generated at the bubble surface is advected and diffused downstream, creating 
a wake whose symmetry is broken by the wall, resulting in a force pushing the bubble away from the 
wall. On the other hand, one expects that a contribution to the lift related to the Bernouilli effect would 
attract the bubble towards the wall. The question is which mechanism dominates. For clean bubbles, 
which generate less vorticity at the bubble-liquid interface than contaminated bubbles, the sign of the 
wall-induced lift changes at Re ~ 35. The bubble shape is known to play a more important role than 
the Reynolds number for the loss of stability of the paths of rising bubbles [37] . The greater the aspect 
ratio between the axis which is perpendicular to the bubble's trajectory and the axis parallel to it, the 
sooner the instability arises, and this is suggested to be related to an increased vorticity generation on the 
bubble surface (see [37]). While this result concerns unsteady flows, we will nevertheless use this finding 
on the bubble shape, as well as the other experimental insights reported in this paragraph, to guide us 
in our own numerical work in Section |4] 

We now introduce the basic mathematical notions for the current work. The Navier-Stokes equations 
in the time dependent domain £l + \B(t), with Q + =Rx[0, oo) and B(t) = {x e fl + \ (x+vgt) 2 + (y—d) 2 < 



r 2 } are 



d t U + U ■ WU + VP - vAU = , 
V ■ U = , 



(3) 
(4) 



3 



where, with x — (x,y) T , U — U(x,t) is the velocity field and P = P(x,t) the pressure field. The 
boundary conditions at the wall (placed at y = 0) and at infinity are 

U\ y=0 = , (5) 
lim [7 = 0, (6) 

x— J-oo 

whereas on the body we may consider either noslip boundary conditions 

U\ gB =v B = (-v B ,0) T , 

or slip boundary conditions 

U ■ n\ dB = , 
t ■ [-P1 + vV{u)\ ■ n\ dB = , 

where n and t are respectively the normal and tangential unit vector at the surface of the body, I is the 
identity matrix and T>(u) — (Vix) + (Vit) T . 

We are interested in solutions which are stationary when viewed in a reference frame comoving with 
the body. In terms of the velocity field u relative to the body, we have 

U{x,t) = u(x - v B t) + v B , (7) 

and u satisfies, in the time-independent domain = fl + \ B, with B = {x E £1 \ x 2 + (y — d) 2 < r 2 }, the 
time-independent equations 

u ■ Vu + Vp - vAu = , (8) 
V • u = , (9) 

with boundary conditions on the wall and at infinity 

U \y = Q = «oo , (10) 

lim U = Uryo , (11) 
\x\— >oo 

with = — v b = { v Bi 0) T - I n terms of u, the noslip boundary conditions on the body become 

u\ aB = 0, (12) 

and the slip boundary conditions become 

u-n\ dB = 0, (13) 
t ■ [-pi + vV{u)\ ■ n\ m = . (14) 



For the numerical treatment of these equations, we solve ([8| and ([9| in the bounded domain fl = Q + \B, 
where (l + = {(x,y) € (— 1,1) x (0,1)} and I > d + r is an arbitrary truncation length. The choice to 
truncate the domain at equal lengths upstream and downstream is motivated by technical reasons, see[B| 
Other choices of domain can be considered, but we do not want to indulge here in questions of domain 
optimization. The truncation introduces artificial boundaries at x — ±1 and y — I. The main focus here 
will be on the choice of boundary conditions on these boundaries, which will be discussed in Section E] 
For convenience later on, we note that the drag and the lift on the body are given by F — (Fd,Fl) , 
where 



F = / (-pI + V(u))ds . (15) 

JOB 

The rest of this paper is organized as follows: in Section[2]we present the different boundary conditions 
which we will implement on the artificial boundaries, in Section [3] we discuss the numerical aspects of the 
work and in Section [4] we present results which validate the adaptive boundary conditions. In Section [5] 
we show that the theoretical behavior of the flow described in |23| is numerically verified in the range of 
simulations we have run. In Section [6] we present the hydrodynamic forces as a function of the body- wall 
distance for different sizes of the body. The appendix, finally, contains technical points concerning the 
adaptive boundary conditions. 
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2 Artificial boundaries 



Theoretically, the correct way to treat the edges of a domain Cl which is obtained by a truncation of the 
half-plane would be to use the solution of the original problem in the half-plane evaluated along those 
edges as a Dirichlet boundary condition. Of course, the solution of the original problem is unknown. 
One must therefore find boundary conditions which represent a good approximation to the solution of 
the original problem. We shall define and investigate three choices: simple boundary conditions (s.b.c), 
classic (or open) boundary conditions (c.b.c.) and adaptive boundary conditions (a.b.c), i.e., the new 
scheme which we propose here. More precisely: 



The s.b.c. simply prescribe on all the artificial edges. 
• The c.b.c. prescribe on the upstream vertical boundary in order to fix the inflow, and impose 



( 14 1 on the remaining artificial boundaries, allowing in- and outflow. 



• The a.b.c. use expressions ( 16 ) and ( 17 1 which are based on the asymptotic expansion of the solution 



of the original problem, to prescribe Dirichlet boundary conditions. 

The s.b.c, while a reasonable starting point since they are in particular also used in the construction 
of weak solutions, are nevertheless problematic, as they do not allow fluid to move through the artificial 
boundary parallel to the wall, making the problem effectively a channel flow. This impacts flow rate 
conservation in two problematic ways: first, the velocity of the fluid must increase artificially above and 
below the body, it cannot be adjusted thanks to fluid "exiting" through the boundary parallel to the wall; 
second, the flow rate should in fact be lower in the truncated domain in comparison with the flow without 
a body, however the use of Uoo at the inlet boundary (i.e., the upstream artificial boundary) prescribes 
the same flow as without a body. See (36| for an example where such boundary conditions are used 
in the three-dimensional version of the problem considered here. A recent work using these boundary 
conditions, albeit for a flow in the full plane around two side-by-side cylinders, is |34| . where the authors 
have run simulations in domains with sizes 750 by 500 cylinder radii to ensure that perturbations due to 
the boundary conditions are small enough. 

The c.b.c. are mixed Dirichlet (pressure) and Neumann (velocity) boundary conditions, and are a 
standard feature of the COMSOL program. They have been used in problems with artificial boundaries, 
see for example [TT] and [53] for the case of an outflow of a channel with a backward-facing step, or again 
as mentioned before, [TJ, for a three-dimensional implementation of the problem considered here. While 
the c.b.c. are less restrictive on the flow rate than the s.b.c, the inlet boundary conditions still prescribe 
a flow rate which is too large. 

We now present our adaptive boundary conditions. As already mentioned before, they are Dirichlet 
boundary conditions on the velocity, based on the asymptotic expansion for the solution of the problem. 
Our adaptive boundary conditions are given by it* = (u* , u* ), with 

«•(*, v) = -oo (i + j-^nWv) + j^^Ax/y) + j^^/y) 

-^*& x/ rt-^*& x/ rt)> (16) 

+ ^ l(4 ^ 2) + ^ 2(4 '^ 2) )' (17) 
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where 



ipx (z 



1 r 

4\/7T 



1 — z 2 + zr + 2z 
r 3 ^/r + 1 



1 r + 1 



2z 



4^ 
1 2z 

1 1 - z 2 
2n r 4 
11-z 2 



r 3 vV + 1 



7r 

1 2z 



(18) 
(19) 
(20) 
(21) 
(22) 
(23) 



where 
and where 



r = yl + z 2 



Vi( z ) = Vw(z) , 
uji(z) = uj w (z) , 

m( z ) = Vb{z) - 2rjw(z) - 2zi]' w (z) , 
lo 2 {z) = cj b (z) - 3uj w (z) - 2ztj' w (z) 



where 



in (26) and (27 1 denotes the derivative, and 

-l/4z 




z>0 
z < 



1/42 



z > 
z < 



2z) 



8nz 



) (e- 1 / 42 -D(l/2|z| 1 / 2 ) 
4z) + v^rH(l - 62) (e- 1 / 42 - D (1/2|^| x / 2 )) ) 



(24) 
(25) 
(26) 
(27) 



(28) 
(29) 

(30) 
(31) 



where 



D(*) = 

is the Dawson function, erf is the Gauss error function 

erf z = —= 








and erfi(z) = — ierf(ijs) is the imaginary error function. 

For a derivation of the adaptive boundary conditions, see [X] In particular, note that we set c*^ = 
throughout the present work. Therefore, the adaptive boundary conditions only reference the charac- 
teristics of the body (size, position, shape, boundary conditions at the interface) through the constant 
c* G M. This constant should be equal to a constant based on the solution in the original unbounded 
domain, defined in (2] and represented there by c\. A good approximation can be determined as part of 
the solution process, for example by the algorithm described in the next paragraph. 

In essence, our algorithm used to determine c\ searches for the root of the function 



9{x) 



- [ V- (T(u,p)V)du>-x 
ni Jn 



(32) 
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where 

T(u,p) = ~u<E)u + vV(u) - pi 

is calculated with u and p computed from the numerical solution obtained with the a.b.c. used with c\ 
such that g(c*) = and where ® represents the dyadic product (i.e., (a (g> b)ij = aibj), and 

^ = (yyxo(^y)>o) T , 

with XCi{ x i y) — Xx(x) ■ Xy(y) a user-defined differentiable cut-off function that cuts a channel perpendic- 
ular to the wall centered on the body and two strips, one adjacent to the wall and one adjacent to the 
artificial boundary parallel to the wall. The scalar n\ in (|32|), given by 



». = / x v m{^T Mz) - Mz) +^ v ^) dz 

Wl *» l V<W/ *> ^ (^1,8)1/4 + £ « (£ 3 /3z) l/4 J ' 

is independent of the solution. Note that, c* and ni depend on I and that for I — ¥ oo the function g(x) 
vanishes for x = c\ {i.e., the exact constant associated to the solution of the original problem), see [B] for 
more details. 

The root-finding algorithm, in our case Brent's method (see [B]), operates on a sequence of simulations, 
refining c\ at each step. Since Brent's method starts with the bisection method, we must first bracket the 
root. We do this by imposing the heuristically derived values '* = 10Re(z — l)d with i = 1,2, . . ., so 
that for the first run the simulation coincides with the one for the s.b.c. As soon as the root is bracketed, 
we begin Brent's method until the desired tolerance is achieved. 



We finally discuss the regime of applicability of our adaptive boundary conditions with regard to the 
size I of the computational domain. There are three conditions that have to be respected. First, Z m i n has 
to be large enough in order for the asymptotic expansion to be indeed evaluated in its domain of validity, 
i.e., y should be large enough for the ratio £ v x/y 2 to be small. Thus, when x ~ y (see [A] for a motivation 
why x large is also a valid region), we must have / m ; n = y ~ C\l v . Second, the artificial boundary 
must be beyond any standing eddies which are not taken into account by our asymptotic expansion. On 
the basis of measurements by Gerrard (see |121 p. 362]), standing eddies typically extend up to three 
times the body diameter at Re « 40, so that we choose ^ m j n = x ~ C^r. Thirdly, the wake should have 
interacted with the wall by the time it reaches the artificial boundary, because the asymptotic expansion 
is made under that assumption. Because the wake scales as £ v x/y 2 and the body is at y = d, we must 
have l min = x ~ C 3 £~ 1 d 2 . Collecting the requirements, one gets 

U ~ maxjd £ v , C 2 r, C 3 £~ x d 2 } . (33) 

Of course, these requirements are qualitative only, since we have no way to estimate the constants C\, C'2 
and C3. They do however explain the observation that the adaptive boundary conditions are particularly 
useful in an intermediate range of parameters, to be specified in what follows. 



3 Numerics 



To numerically solve (J8JH 10 ) with all the different choices of boundary conditions, we use COMSOL 
Multiphysics 3.5a, controlled through a Matlab 2009a script. The linear system is solved using a direct 
linear solver, PARDISO, and the nonlinear solver is a damped Newton fixed point solver. The mesh 
elements are a mix of triangles and quadrilaterals of Lagrange type, of degree two for the velocity and 
degree one for the pressure. For all the simulations presented below, a workstation with 36 GB RAM and 
a 12-core processor was used. The algorithms mentioned above are features of the COMSOL program, 
and we refer to the product documentation for additional information. 
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4 Validation 



The first set of numerical tests is used to validate the adaptive boundary conditions. For the sake of 
convenience, we have set 2r = 1 for all simulations of this validation run, so that the Reynolds number is 
given by the inverse of the viscous length. The body is at a fixed distance d = 1 from the wall. We show 
the effect of the domain size for all choices of boundary conditions on the computation of the drag and 
lift. 

4.1 Mesh 

We now present our method for generating meshes for the domain Cl. Figure [l] represents the coarsest 
and smallest mesh used, from which all others are constructed. The smallest mesh will have a domain 
with truncation length I = 10 in the same units as the body size. 



Figure 1: The coarsest mesh (h = 1.25) on the smallest domain (/ = 10) used in the simulations of 
Section |U 

We adopt the same philosophy as [17 : in order to minimize possible effects from one mesh to another 
due to element placement, we define a small rectangular zone (of size 10 by 5) around the body meshed 
with a constant node distance h on the top, left and right boundaries, and using triangles placed according 
to an advancing front algorithm provided by COMSOL. Since the body size is fixed, this mesh will be the 
same for all simulations. The rest of the domain is paved with strictly identical square elements with sides 
of length h. Increasing the domain size then consists in adding supplementary square elements of the 
same size, guaranteeing that the meshes always have the same structure regardless of domain size. The 
base mesh for I = 10 has ho — 1.25, 367 triangles and 96 squares, for a total of 2718 degrees of freedom. 
A step of mesh refinement simply consists in dividing all mesh cell lengths by a factor m r — 2 effectively 
dividing each element into four identical smaller ones. Thus the recurrence relation hi — hi—i/m r . 

We compute simulations with mesh sizes h to h 3 on domains with truncation lengths I = 10, 20, . . . , 90. 
All figures in this article representing the velocity fields are obtained from simulations computed on meshes 
with elements of size /13. Since the drag and the lift are given by integrals, we take advantage of the 
fact that each of our successively refined meshes is simply a subdivision of the preceding one and use 
a Richardson extrapolation scheme (see [28 ) to accelerate convergence. We base our scheme on the 



hypothesis that the error of the drag and the lift (151 is of the form 



h n =I+ C 2 h 2 n + C 3 h 3 n + C A h 4 n + 0(h 5 n ) , (34) 
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with Ci unknown constants. This yields the extrapolation formula 



ml ■ I hs - (ml + + m 5 r ) ■ I, l2 + (mf. + mf. + m 2 ) ■ I h - I h() n(h5 s 

+ U(n ) 



(ml + m% + ml) + (mj. + mf + m 2 ) — 1 

0(ht) , (35) 



512 • I h3 - 224 ■ I h2 + 28 ■ I hl - I ho , 5 



315 



for the Richardson extrapolate Jr of /. The error terms in ( 34 ) start at because the Lagrange elements 
are of degree two for the velocity and one for the pressure. Note that if our Ansatz (34) should miss 
intermediary terms of order h n , n > 2, then they are not eliminated by the scheme, but only somewhat 
damped. Thus, the extrapolation scheme will always diminish the magnitude of the error, even though 
not necessarily to order 

4.2 Results 

4.2.1 Qualitative aspects 

In a first step, we show that the choice of boundary conditions has an important impact on the quantities 
being calculated. To showcase this, we represent in Figure [2] the streamlines for the velocity field as 
seen from the body, i.e., the velocity field from which the constant flow has been subtracted. In 
the case of the simple boundary conditions (s.b.c), we observe an artificial backflow (moving clockwise 
on our figure) . A very important proportion of the computational domain is thus used to compute this 
non-physical flow. In the case of the classic boundary conditions (c.b.c), this backflow is not present, 
but the flow is still significantly influenced by the artificial boundaries, as can be seen by looking at 
the streamlines. The adaptive boundary conditions (a.b.c.) are the only ones that yield qualitatively 
satisfying flows all the way up to the boundary, and the streamlines exhibit almost no distortion near the 
boundary. Figure [3] presents the streamlines for Reynolds numbers Re = 0.5, 5, 10, 25. 

We now discuss the qualitative behavior of the flow computed with the a.b.c. in more detail. The 
case of Re = 25 shows signs of artificial behavior at the right-hand boundary, which could be due to 
the fact that the asymptotic expansion ignores the actual position of the body. More precisely, the 
a.b.c. assume that the wake has already interacted with the boundary before the position of the artificial 
boundary. The wake is asymptotically governed by the first order in the £„a;/?/ 2 -scaling of the asymptotic 
expansion, and larger Reynolds numbers therefore yield narrower wakes which interact with the wall 



farther downstream, see (33 1. In order to simulate higher Re flow, one would have to bring the body 
closer to the wall, which would change the actual problem being solved, or use a larger domain, which 
was outside our computational power. This limitation to Re < 25 is not so stringent, because in the 
absence of any wall, the flow behind a circular body becomes unsteady at Re ~ 30, see |35[ p. 150], and 
we do not expect that the situation be radically different in our problem. 



The case of Re = 0.5, also shows distortions in the streamlines, especially downstream. Since Re = £ v 



-l 



in the validation run, flows with Reynolds numbers smaller than unity need a large domain too, see (33 1. 
Indeed, since the wake scales like £ v x/y 2 , a low Reynolds number actually delays the distance after 
which the solution may be represented by the asymptotic expansion (which is obtained for y — > oo, i.e., 
£ v x/y 2 — > 0), so that the a.b.c. should be used in conjunction with larger numerical domains. As can be 
seen in Figure [4] increasing the domain size does indeed improve the behavior of the streamlines. This 
case is different from the large Reynolds number case, since the wake does interact with the wall, but the 
flow region which is represented is too small compared to the viscous length for the asymptotic expansion 
to be accurate (the flow would probably be more accurately represented under the Stokes approximation) . 
It is nevertheless possible to simulate flows for lower Reynolds numbers, but one must respect particular 
scaling conditions. This procedure will be presented in Section [5] where we will investigate a sequence of 
simulations with progressively smaller bodies. 

4.2.2 Quantitative analysis 



We present in Figures [5j|9] the drag and lift for a sequence of simulations with Reynolds numbers Re = 
0.5, 1, 5, 10, 25, performed with the three boundary conditions. As the simulations with adaptive boundary 
conditions overall present the least variation with domain size, we use the largest simulation (I = 90) 
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Figure 2: Streamlines for a domain with truncation length I = 10 and Reynolds number Re = 1, for a 
body with noslip boundary conditions. The colored map represents the velocity norm. Streamlines (top 
left). Streamlines as seen from the body: s.b.c. (top right), c.b.c. (bottom left) and a.b.c. (bottom right). 
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Figure 4: At Re = 0.5, a larger domain diminishes the distortions in the streamlines. 



feasible on our workstation using these boundary conditions as a reference for the computation of relative 
errors. 




Figure 5: Re = 0.5 simulations, for a body with noslip boundary conditions. Top: Drag (left) and lift 
(right) as a function of domain size for the three boundary conditions. Bottom: relative error on drag 
(left) and lift (right) as a function of domain size (log-log scales). 



In Figures flQ^fTT] we summarize the relative errors for the drag and lift obtained for a body with slip 
boundary conditions. 

In Tables [T]-[2] we present the extrapolated values of the drag and lift for circular objects with noslip 
and slip boundary conditions at Reynolds numbers Re = 0.5,1,5,10,25. The extrapolation was done 
using formula ( 35 1 on the values obtained from simulations with mesh sizes ho to I13 , for the largest 
domain (Z = 90). 

We observe that, for both boundary conditions on the body, the simple boundary conditions overes- 
timate the values of the forces compared to the classic boundary conditions which in turn overestimate 
the drag and lift compared to adaptive boundary conditions, except for the drag in the case Re = 0.5 
on the smallest domain. This was to be expected, since both the s.b.c. and c.b.c. impose a flow rate 
only appropriate in absence of a body. In a bounded domain, this flow rate is higher than when there is 
a body, leading to an inevitable overestimation in the hydrodynamic forces. However, the modification 
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•Noslip body' 





Re = 0.5 


Re = 1 


Re = 5 


Re = 10 


Re = 25 




s.b.c. 


20.020 


10.601 


2.9595 


1.9164 


1.1896 


drag: 


c.b.c. 


20.011 


10.597 


2.9582 


1.9155 


1.1888 




a.b.c. 


20.007 


10.594 


2.9571 


1.9145 


1.1880 




s.b.c. 


1.4480 


1.3912 


1.0454 


0.77384 


0.37497 


lift: 


c.b.c. 


1.4460 


1.3898 


1.0448 


0.77344 


0.37481 




a.b.c. 


1.4450 


1.3890 


1.0442 


0.77304 


0.37464 



Table 1: Extrapolated values for drag and lift computed with the largest domain (7 = 90) for a body 
with noslip boundary conditions. 



"Slip body" 





Re = 0.5 


Re = 1 


Rc = 5 


Re = 10 


Re = 25 




s.b.c. 


14.505 


7.6636 


2.0929 


1.3023 


0.71367 


drag: 


c.b.c. 


14.500 


7.6611 


2.0922 


1.3018 


0.71338 




a.b.c. 


14.497 


7.6597 


2.0916 


1.3014 


0.71307 




s.b.c. 


0.88168 


0.84854 


0.59285 


0.37608 


0.10188 


lift: 


c.b.c. 


0.88065 


0.84779 


0.59250 


0.37585 


0.10176 




a.b.c. 


0.88013 


0.84732 


0.59219 


0.37563 


0.10162 



Table 2: Extrapolated values for drag and lift computed with the largest domain (I = 90) for a body 
with slip boundary conditions. 
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Figure 6: Re = 1 simulations, for a body with noslip boundary conditions. Top: Drag (left) and lift 
(right) as a function of domain size for the three boundary conditions. Bottom: relative error on drag 
(left) and lift (right) as a function of domain size (log-log scales). 
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Figure 7: Re = 5 simulations, for a body with noslip boundary conditions. Top: Drag (left) and lift 
(right) as a function of domain size for the three boundary conditions. Bottom: relative error on drag 
(left) and lift (right) as a function of domain size (log-log scales). 
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domain size domain size 



Figure 8: Re = 10 simulations, for a body with noslip boundary conditions. Top: Drag (left) and lift 
(right) as a function of domain size for the three boundary conditions. Bottom: relative error on drag 
(left) and lift (right) as a function of domain size (log-log scales). 
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Figure 9: Re = 25 simulations, for a body with noslip boundary conditions. Top: Drag (left) and lift 
(right) as a function of domain size for the three boundary conditions. Bottom: relative error on drag 
(left) and lift (right) as a function of domain size (log-log scales). 
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of the flow rate due to the body decays as the inverse of the square root of the size of the domain, so 
that in the limit of a complete half-plane, the difference vanishes. This is different from the case of a 
body in the whole space, where the presence of abody imposes a non-vanishing modification to the flow 
rate (see [I] for the asymptotic expansion for the velocity field in this case). Nevertheless, in the case of 
truncated domains we show that the difference in flow rate implied by the various boundary conditions 
has a non-negligible effect. 




- classic 

- simple 

- adaptive 



30 40 
domain size 



30 40 
domain size 



50 60 70 
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adaptive 
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Figure 10: Simulations for a body with slip boundary conditions. Relative error for the drag as a function 
of domain size (log-log scales), for Re = 0.5, 1, 10 and 25. 



4.3 Body shape 

We next demonstrate the range of applicability of the a.b.c. All simulations were performed using a triple 
refinement of the mesh shown in Figure |12| obtained according to the advancing front algorithm provided 
by COMSOL (maximum triangle side is 1.0). We recall that one step of refinement is defined as the 
division of all triangles into four smaller ones of equal size. The characteristic body size is defined by 
the projection of the shape on the y-axis. We present four shapes in Figure [T3| an ellipse with its large 
axis perpendicular to the flow of length 1.0 and height 2.0, a square of side a/9/8, a "marmite" composed 
of half a circle and half an ellipse whose longer axis is twice the smaller one, of total height 1.0, and a 
"bone" of length 1 + \[2 and height 1.0. These bodies serve to illustrate that the a.b.c. work well for 
shapes with concavities and asymmetries, as well as non-smooth boundaries, even though the theoretical 
developments are expressly derived for smooth bodies. 

It is also possible to simulate collections of bodies. We tested arrangements of two and three circles, 
contained in a circular area of radius r — 1.0, see Figure |14| These two simulations were done on a 
smaller domain, with I = 10 (the mesh was generated in the same way as for the other shapes). 

5 Conformance to expected theoretical behavior 

A result of |2"31 Section 2.3] states that solutions of (JsJ) and ^ tend to zero when the obstacle size tends to 
zero. This theoretical result may seem obvious at first, but because of the Stokes paradox this statement 
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domain size domain size 



Re = 10 Re = 25 




domain size domain size 



Figure 11: Simulations for a body with slip boundary conditions. Relative error for the lift as a function 
of domain size (log-log scales), for Re = 0.5, 1, 10 and 25. 



Figure 12: Mesh composed of triangles obtained from COMSOL's advancing front algorithm on a domain 
with truncation length Z = 20, a circular body with r = 0.5 and d — 1.0. 



requires a proof. Since the FEM resolution scheme used by COMSOL is based on a Galerkin method, 
it is natural to test this result. In the present work we examine the relative velocity, so that the result 
shows that for vanishing body size, the flow tends to u^. 

Having validated our scheme, we exclusively use the adaptive boundary conditionsfrom now on. In 
what follows, the viscous length is fixed for all simulations at £ v = 1, all the bodies have noslip boundary 
conditions, and their sizes are chosen as 2r = 0.005, 0.01, 0.05, 0.1, 0.5, 1, 2. For reasons of computational 
tractability and because of the requirement of placing the artificial boundary far enough to ensure that 
they are evaluated in the asymptotic physical regime, we choose the truncation length of the domain to 
be I — 20 for all simulations. This size is already quite impressive compared to the size of the bodies. 
Since we work with one fixed domain size only, the meshes are for simplicity chosen to be triangular 
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Figure 13: Top left: ellipse; top right: square; bottom left: "marmite"; bottom right: "bone". 
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Figure 14: Collections of two and three circles. 

pavements of the fluid domains, as given by COMSOL and similar to that of Figure [12] We however go 
as far as quadruple refinement (leading to approximately 3 million d.o.f.) and then apply an appropriate 
Richardson extrapolation, similar to ( |35[ ). Meshing difficulties arise with smaller bodies, in the sense that 
the mesh elements become smaller near the body, thus posing a problem with regard to memory or mesh 
quality (ratio of the smallest to the largest element). 

Following the proof in |23| , we wish to provide evidence that the norm of the velocity gradient inte- 
grated over the fluid domain vanishes as the body disappears. The norm of the velocity gradient is given 
by 

||Vu|| = yj(d x u) 2 + (d yU y + {d xV y + {d y vf . (36) 

For the sake of comparison, we also compute the norm of the velocity field as seen by the body, i.e., 

\\U - UooW = y/(u - Uoo) 2 + V 2 , (37) 
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and both norms are integrated over the whole computational domain and normalized by its surface. 



Figure 15 shows the evolution of these norms as a function of body diameter (2r). Over the range 



of body sizes that we investigated, the results suggest that the norm of the velocity gradient vanishes as 
| log as r — > 0. 



10 - 




10 ' 



10' 

Body diameter (2r) 
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Figure 15: The norms of the velocity gradient and the velocity field as seen by the body, as a function of 
body size normalized by the fluid domain surface. 

This investigation shows that given the right choice of parameters, the adaptive boundary conditions 
may be used to simulate flows with smaller Reynolds number (in this run Re = 2r) than what was 
expected at first from the results of the validation. Indeed, here the body size is decreased and the 
viscous length remains the same (resulting in a decrease of the Reynolds number), and thus the a.b.c, 
which explicitly depend on £ v , are always evaluated at some well-adapted distance regardless of body 
size, whose effect is only reflected through the value of c\. This is in contrast to the simulations presented 
in Section [4j where the body size was kept fixed, and the Reynolds number (and thus the viscous length) 
changed, which would have required to adapt the domain size appropriately for the a.b.c. to be useful. 

6 Force as a function of wall distance 
6.1 Mesh 

In this investigation, we keep to the triangle-only family of meshes, similar to the ones already used in 
Sections |4 . 3 1 and |S"| We choose a domain truncation length I = 40 for all simulations, and the most refined 
meshes are composed of almost 370'000 elements and approaching 1.7 million degrees of freedom. 



6.2 Results 

We only treat the case of bodies with slip boundary condition, because, according to [30 , it is this type 
of body which is susceptible to experience a zero transverse force for a certain body- wall distance. We 
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present the results for simulations for bodies with circular and elliptical shapes. For the case of the 
elliptic bodies, we choose a low ratio of axes a x /a y , since this is the shape with the smallest lift for a 



given body- wall distance (see Figure 16 1. An ellipse with its longer axis perpendicular to the flow direct 
may also act as a simple model for bubbles in flows of Reynolds number Re ~ 0(100), where real bubbles 
flatten, see |57] , 

Again we tested the range of Reynolds numbers Re = 0.5, 1, 5, 10, 25 (obtained numerically by changing 
the dynamic viscosity v), for the range d=0.6,...,2.5 (with increments in steps of 0.1). Figure 17 shows 
the hydrodynamic forces in the case of a circular body at Re = 1, and Figure [18] shows the same forces 
for the elliptic body with a x /a y — 0.1 at Re = 25. In none of these combinations of parameters do we 
observe a change of sign in the lift, in contrast to what is predicted by the experiments for the three 
dimensional case. 




0.8 1 1.2 

body aspect ratio (ax/ay) 



Figure 16: Hydrodynamic forces at Re = 1 versus ellipse aspect ratio, for a body to wall distance d = 1. 
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Drag force in function of body position, slip bubble. Re = 1 
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Lift force in function of body position, slip bubble. Re = 1 
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Figure 17: Hydrodynamic forces versus body-wall distance, for a circular object, at Re = 1. 
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Drag force in function of body position, slip bubble. Re = 25 

i i i i i i i i i i i r 
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Figure 18: Hydrodynamic forces versus body-wall distance, for an elliptical object with aspect ratio 
a x /ciy = 0.1, at Re = 25. 

7 Conclusion 

We have presented simulations of flow around bodies in a half-plane for a wide range of parameters, 
including a range of three orders of magnitude of Reynolds numbers for some cases. We have shown 
that thanks to the adaptive boundary conditions, the dependence of the computed values for drag and 
lift on the size of the computational domain is drastically reduced, achieving accuracy better by one 
to two orders of magnitude when compared to simulations with simple or classic boundary conditions. 
Therefore, with the a.b.c. a given accuracy can be obtained on much smaller domains, thus bringing down 
the hardware requirements (CFD on a laptop). We have also shown a substantial qualitative improvement 
of the physical behavior of the flow, in the sense that the adaptive boundary conditions have a minimal 
influence on the streamlines, in particular close to the artificial boundaries. 

The results obtained in the numerical simulations seem to indicate that the lift monotonically decreases 
with the distance to the wall, but no point where the transverse force vanishes could be found. It is an 
open question if this experimental result depends on dimensionality, on the exact geometry of the body 
or fluid container, or if it is the result of the flow being unsteady. 

Finally, the adaptive boundary conditions has been shown to work well for a variety of smooth bodies, 
even in collections. 
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A Deriving the asymptotic boundary condition 



We show how to derive the adaptive boundary conditions ( [161 ) anc ^ (17 1. In [23 it was shown that the 
non-dimensional system 



dxii + u ■ Vtt + Vp — Am = , 
V • u = , 



(38) 
(39) 



in the domain Q, = {(x,y) € K x [l,oo) \ £>}, is equivalent to the same problem without the body, but 
with some force term of compact support, in the sense that the solutions coincide outside some compact 
set containing the body. For the system without a body, an asymptotic expansion of the velocity field 
(with y^ 1 playing the role of the small parameter) was obtained in [2], given by 



v lls (x,y) 



^cpx(x/y) + ~<p2,i(x/y) + 4^2,2 (a/ jr) 

yJ/Z yZ yZ 



^Vw(x/y ) - ^r) B {x/y ) , 
1 ipi(x/y) + ^ip2,i{x/y) + ^2 t i{x/y) 



f.3/2 



+ —u w (x/y ) + —uj B (x/y ) , 

yi yi 



(40) 



(41) 



where (x,y) elx [l,oo), and where the functions ipi, op2,i, P2.2, ipi, "02,1, "02,2, Vb, Vw, and low are 
as defined in (Jl8|-(23l and ( 28 )— ( 31 1 . The constants c\ and C2 depend on the solution and are defined 
more precisely in [2] . The equations ^ and (J9j) are obtained from ( 38 1 and ( 39 1 by setting 



u(x) = u(x) — ei 



x 

C^ + i 



(42) 
(43) 



Inserting this into ( 40 1 and (|41| we get 



u < J-u aa (x,y) = 1 + u^ix/iv, 1 + y/i v )) 

Cl 



1 + 



For y — » 00 we have, for example 



and 



r)W 



x I ' £ v 



(1 + y/4) 3 / 2 

Cl 



<P2,i{x/y) 



(i + y/L) 2 y 



2 



y3 



C"2 



(i + y/4) 2 



(i + y/4) 2 

Cl 



V>2,2{x/y) 
VB(x/y 2 ) 



Vw 

Vw 
Vw 



( x/l v 

\(y/Q 2 



2-^?^ X / 



y y 2 



-Vw 



(44) 



where we have used the Taylor series of the function rjw for small £ v x/y 2 . Applying this to each term 
in (44 1, then discarding terms decaying faster than y~ 2 obtained from the x/y scale, and decaying faster 
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than y 3 obtained from the x/y 2 scale, and remembering that any explicit x be grouped with scaling- 
appropriate powers of y~ 1 , we get 

«M«as(Z) y) ~ 1 + u as (x/£ v ,y/£ v ) - 2 ^ L 'n w (i v x/y 2 ) - ^^^^^-r/^^x/y 2 ) < 

and similarly, this time discarding terms decaying faster than y~ 4 obtained from the x/y 2 scale 



Uoo^asfoj/) ~ v as {x/l v ,y/l v ) 



uj w (£vx/y 2 ) - 



uj w (l v x/y 



T v y 

This procedure is a resummation technique, which happens to yield a new expansion which respects the 
boundary condition (10 1 if we set ci — 0, i.e., 



lim u as {x, y) = u(x ^ 0, 0) = m c 



Remark 1 is </ie fact that the asymptotic expansions involve scale-invariant functions which makes 
this exchange of limits possible, allowing to extend the validity of the expansion, which was originally 
valid for large y only, also to large x. A detailed mathematical proof of this feature has however not been 

) and now follows from these results, with the constants c\ and d^, 
replacing c\ and C2 since the domain is finite in the numerical implementation (see\B\for more details). 



carried out yet. Expressions (16 



B Approximating the unknown constant of the asymptotic ex- 
pansion 



We motivate the choice of the function g given in (32 1. First, we must find an asymptotic description for 



the pressure, valid in the domain fl. Using the Ansatz 

1 

P = ^2 Ue ' Ue + 9 ' 

where the index E means that we only retain functions with x/y as the argument in the asymptotic 
expansion of the velocity. Inserting this into ^ we get 

u ■ Vw — ue ■ Vme; + Vp — isAu = , 

and then taking the divergence yields 

Ap = - (Vu) T : (Vu) + (Vm b ) T : (Vu B ) , 

where " : " denotes the tensorial scalar product (i.e., A : B = £^ • a#6y). The dominant terms on the 
r.h.s. are products of functions of x/y and x/y 2 , or x/y 2 and x/y 2 , decaying at least as fast as y~ lx l 2 . 
Since for large y the dominant terms on the l.h.s. are those obtained from d 2 , integration shows that 
p y~'' 7 1 2 . which is beyond the terms we retain, so that to leading order the pressure is given by 

p ~ - ul G + ww^ Mx/y) + wh? V2Ax,v) ) ■ (45) 

Second, we obtain a formula which gives the constant c\ as a ratio between an integral over the domain 
f] and an integral over its boundaries. We define the tensor T by 

T = -u <g) u + vV(u) - pi , 

where ® represents the dyadic product (i.e., (a ® b)ij — aibj) and T>(u) = Vit + (Vu) T , with V • T 
yielding We choose the vector 

v = (Vyxn(x,y),o) T , 
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where the factor y/y is used for reasons that will become clear later on and where 

Xn(x,y) = Xx(x) ■ xy(y) 

is a cutoff function which cuts a channel perpendicular to the wall centered on the body of radius r and 
two strips, one adjacent to the wall and one adjacent to the artificial boundary parallel to the wall of the 
domain f2 = Cl + \ B, where Cl = (—1, I) x (0, I) and B = {x e Cl \ x 2 + (y — d) 2 < r 2 } (d is the distance 
between the body center and the wall), with I > d + r and d > r > 0. Namely, 

Xx(x) = X (-{r + x)/A) ■ X ((x - r)/4) 

and 

x y (y) = x(y/4)-x((i-y)/4) , 

where \ i s an arbitrary smooth cutoff function which we have chosen as follows, 



^ /■max{0,min{l,?7}} 

X(V) = - V(l-V) 4 d(, 
P Jo 



dvXiv) 



, otherwise 



and p = Jq r](l — r/^dr]. The factor 1/4 in the arguments of \x and \y IS arbitrary and is chosen such as 
to have numerically reasonable gradients. The power 4 in the definition of x IS arbitrary as well and is 
simply chosen to yield a sufficiently smooth expression. We have 

Tv=v»»(^)(-(::) + ,( ai x)-(o))- 

Integrating over the whole domain, we get, by Gauss's theorem, 

/ V • TVdu = [ (TV) ■ nda + [ (TV) ■ nda . (46) 
Jn Job Jon 

The integral over the surface of the body dB vanishes by the choice of the cutoff function. We have 
:= I V • TVdu 



Xx(x) \ ^ V ^ + VydyXyiy)) (-uv + vd x v + ud y u)duj 
+ / VyXy(y)dxXx(x)(-u 2 + 2vd x u-p)du: , (47) 



which is computed numerically from the FEM solution. Due to the choice of the cutoff function, we also 
have 



hn I (TV) • nda = [ 
Jon Ji 



VyXn(x,y)(-u 2 +2i/d x u-p)\ (-dy) 



+ 1 ^yxn{x,y)(u 2 -2vd x u+p)\ x=+l dy . (48) 



On dCl we use (16 1 and (17 1 to represent the velocity field. We do not consider terms which decay faster 
than 1/y 2 in the x/y scaling, since these would be of the same order as those we neglect in our asymptotic 
expansion. In the x/y 2 scaling, we neglect those terms which decay faster than 1/y 3 for the same reasons. 
Mixed terms (comprised of a product of functions of either scaling behavior) are neglected if they decay 
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faster than 1/y 2 , meaning that none are actually retained. We obtain 



(y/l v ) 2 ,iy " ,a ' (y/£ v ) 



-2d x u ~ 



P~- \( u e + v e) 



,2 



so that 



(w 2 - 2d x u + p)\ x=+l - (u 2 - 28 x u+p)\ x= _ l 

- ^fm(ivi/y 2 ) + o - ^ {mihi/y 2 ) mi-tJ/y 2 )) 

y y 

- ^miU/y 2 ) - ^ M4*/y 2 ) - m(-tJ/y 2 )) , 

y y 

where we use that ipx(—z) = tpi(z) and that </?2.i is an odd function. We insert this approximation 
into (48 1, where Xai^hv) = Xy(y) («-e., the fact that the computational domain is of equal length 
upstream and downstream simplifies the expressions), and treat the integral separately according to the 
two scalings. We get 

I aii :=c iV / XuU/) dy + 2c x t v J Xy{y)^V2,i{l/y)dy 

^ jf 00 Xjf( i/,) ^ /2 ^M^iM + 2*^) * , 

where we see that the factor ^/y in V is chosen in order to obtain a non- vanishing expression for large Z, 
and 

Note that in the derivation of the integrals on the upstream and downstream boundary of the truncated 
domain Cl, we assumed that the asymptotic expansion may be extended, for large fixed x, to all y although 
the asymptotic expansion is a priori only valid for large y. See Remark [T] in [X] for a motivation. 
We then define 

m(l,e v , XCl ) := (I^+I^/ycl, 

which can be computed numerically. We then have that 

m 

and we expect that the constant c\ associated to the solution of the original problem is given by 

C\ = lim c*(l) . 

I— ^oo 



This motivates the definition of ( 32 ) . 
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